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Abstract 

In this work we investigate the issue of non-physical slip at wall of lattice Boltzmann simulations 
with the bounce-back boundary scheme. By comparing the analytical solution of two lattice models 
with four and nine discrete velocities for the force-driven Poiseuille flow, we are able to reveal the 
exact mechanism causing the issue. In fact, no boundary condition is defined by the bounce-back 
scheme for the the discrete velocities parallel to wall. Other factors, such as initial conditions and 
inlet and outlet boundary conditions, can play the role and induce the non-physical slip velocity. 
Therefore, the issue is not related to the single-relaxation-time scheme. Naturally the key for 
resolving it is to specify the definition for these velocities. Through a lid-driven cavity flow, we 
show that the solution can be as easy as no extra effort required for simple geometries, although 
further study is necessary for complex geometries. 
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I. INTRODUCTION 


The lattice Boltzmann method (LBM) has been developed as a mesocopic computational 
fluid dynamics (CFD) tool for the Navier-Stokes (NS) level problems and beyond jIH3|. Due 
to its origin from the lattice gas automata (LGA) [3], it keeps the flexibility of a particle 
method to a great extent. On the other hand,the stochastic noise is eliminated in LBM by 
using the distribution function. Importantly, this links the LBM into the discrete velocity 
method of the Boltzmann-BGK (Bhatnagar-Gross-Krook) equation[5H8]. From such a point 
of view, we actually solve a system of partial differential equations with linear advection 
terms. This opens the door of introducing more sophisticated scheme leading to such as 
finite difference LBM or finite volume LBM (e.g., [9] and TO- 

It is the simplicity which brings the popularity of LBM. The algorithm is easy to un¬ 
derstand for application purpose. Since an explicit scheme is employed, the programming 
and parallelism is straightforward. The second order accuracy in both space and time is 
achieved at the expanse of a first order scheme, which is fairly enough for most purposes. 
The boundary treatment, even for complex geometry, can be incredibly simple due to the 
so-called bounce-back (BB) scheme which only requires the particles to reverse their velocity 
on the wall/obstacle mm . 

However, there are non-physical slip velocities occurring at wall in simulations using the 
BB scheme. This was firstly discovered for two LGA models HU and then was analysed 
for the nine-discrete-velocity (D2Q9) LBM in [T2]. By using a simple force-driven Poiseuille 
flow, it was shown that the non-physical slip can be generated on the wall with the BB 
family scheme ra- Since the slip velocity was found to be related to the mesh size, it was 
then deemed as a numerical artificial effect. Later, this issue has been considered as an 
inherent deficiency of the single-relaxation-time (SRT) scheme since it may be resolved by 
using extra free parameters in two-relaxation-time (TRT) or multi-relaxation time (MRT) 
schemes (see e.g. USD- Indeed, it is believed that the SRT plus BB combination cannot 
avoid this issue m- 

However, due to its simplicity, the SRT plus BB combination is more favourable for 
application purpose. Therefore, it is of interest to investigate how the slip velocity is induced 
and therefore gain useful information on how to fix it. For this purpose, we will first analysis a 
lattice model with four discrete velocities (D2Q4) for the force-driven Poiseuille flow following 
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the method presented in [12]. With this lattice model, we will see that the SRT plus 
BB combination does not necessarily induce non-physical slip velocity and it is possible to 
correctly implement the non-slip wall. Then we will compare this model with the D2Q9 
model to find the exact mechanism inducing the slip velocity. With these findings, we 
devise a guidance on how to fix the non-physical slip velocity. Finally, we will examine this 
guidance and thereby the discussions on the mechanism by simulating the lid-driven cavity 
flow with the D2Q9 lattice. 


II. LATTICE BOLTZMANN SCHEME AND LATTICES 


The LBM can be considered as an approximation to the Boltzmann-BGK equation [M- 
After the discretisation in the particle velocity space, the governing equation becomes 


df 1 

+ c a • v/ Q = ~-(/ a - r a q ) + f w 


( 1 ) 


which represents the evolution of the distribution function f a (r,t ) for the ath discrete ve¬ 
locity C a at position r = (x, y, z ) and time t. The effect of external body force is described 
by F a . The particle interaction is modelled by a relaxation term towards the discrete equi¬ 
librium distribution function In order to simulate incompressible and isothermal 

flows, it is common to use an equilibrium function with second order velocity terms, i.e., 


fa = w a p[ 1 + 


U ■ C a 1 (U -c ( 


+ — 


RT 0 2 (RT ( 


u-u 

2 RTn 


( 2 ) 


which is determined by the density, p, the fluid velocity, U, and the reference temperature, 
T 0 . For gas flows, the constant, R , can be conveniently understood as the gas constant. If 
a liquid fluid is involved, it, together with To, can be considered as a reference quantity. 
For convenience, the sound speed c s is often considered equal to \ZRTq, although there is a 
constant factor of difference. The weight factor is denoted by w a for the discrete velocity 
C a . The term $F_\alpha$ can be obtained by using various method [HUTS]. Here, the first 
order expansion is sufficient for our purpose, which can be written as [8j. 


F a = pw ( 


G Cg 

? RTn 


( 3 ) 


where the actual induced acceleration is denoted by G. The relaxation time , r, is related 
to the fluid viscosity , /r, and the pressure, p, via the Chapman-Enskog expansion, i.e., 
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/i = pr. Hence, for isothermal and incompressible flows, the Reynolds number becomes 
Re = p 0 U 0 L/h = U 0 L/(tRT 0 ), where we use a subscript 0 to denote the reference value and 
L the characteristic length of the system. It is worth noting here that the Knudsen number 
can be defined as /io \JRTq / (poL). So, the relaxation time r is also related to Knudsen 
number by the viscosity Kn = t\/RTq/L, where po = PqRTq is applied. In this sense, we 
have Kn x Re = Uq/\JRTq = Ma. To get the density and velocity, we only need summation 
operations, i.e., 

P = ^2fa, and, P U = Y. f aCa ' 

a a 

To numerically solve Eq. ([Tj) , a smart trapezoidal scheme can used to achieve the particle- 
jump like simulation |16j, which can be written as 


fair + C a dt, t + dt) - f a {r, t ) 


dt 


t + 0.5 dt 


fa(r,t) - fa q (r,t) 


t dt 

+ ( 4 ) 

r + 0.5 dt 


where 

/« = /« + §(/« - C) - f F a . (5) 

By using f a this scheme is ready for implementing the stream-collision algorithm. At the 
same time, the macroscopic quantities become 


P 


^2 fa, and, pU = 22 C cJ + 


pGdt 

2 


( 6 ) 


For two dimensional flows, the D2Q9 lattice is commonly used where the nine discrete 
velocities (a = 1..9) are 


C a , x = VWTo[ 0,1,0, -1,0,1, -1, -1,1], 


( 7 ) 


C a , y = V3RTo[0, 0,1,0,-1,1,1,-1,-1], 

and the corresponding weights are 


4 11111 1 1 1 

9’ 9’ 9’ 9’ 9’ 36’ 36’ 36’ 36 J 


( 8 ) 

( 9 ) 


As discussed above, the stream-collision algorithm is ready to be implemented now. The 
only trick is to tie the space and time step together as dr = C a dt. For instance, assuming 
the system length is L, we may set the spatial step dx = L/N and then dt = L/(N^/3RT 0 ) 
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Figure 1. Illustration of the D2Q4 and D2Q9 lattice at the bottom wall for a force-driven Poiseuille 
flow. The discrete velocity C\ = (0,0) of the D2Q9 lattice is not shown. 


where N is the cell number. This insures the “particles” are jumping on a uniform grid 
system. In simulations, it is common practice to use a non-dimensional system in which 
the space and time step are considered as reference value. Apparently, this will make no 
difference on results. However, confusion may be caused in this way. We shall return to 
this point below. Alternatively, we may also transform Eq. 0 to its non-dimensional form 
first by using the reference values presented in P2] and then apply the scheme Eq. Q. 
Again, this non-dimensional transformation will not alter the final simulation results but 
the relations with dimensional quantities are more clear, at least for gas dynamics. 

To study the effect of the BB scheme on the solid boundary, we will first use a D2Q4 
model P| where the four discrete velocities are 


C a , x — \J AT 0 [1,1, —1, —1], 


( 10 ) 


C a , y =yfKT Q [- 1,1,-1,1], 


( 11 ) 


and the weights are 


w a 


1111 
^4’4’ 4’4 


( 12 ) 


Both the D2Q9 and D2Q4 models are illustrated in Fig. [I] It is worth noting that in the 
D2Q4 model there is no discrete velocity parallel to the wall for the force-driven Poiseuille 
flow with regular shaped channel (i.e., the wall is either horizontal or vertical), which will 
make a dramatical difference from the D2Q9 model for the force-driven Poiseuille flow. 
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III. SLIP VELOCITY AND BOUNCE-BACK SCHEME 


A. Solution of D2Q4 model for force-driven Poiseuille flow 

In the following, we will use the method presented in ca to solve the D2Q4 model 
for the force-driven Poiseuille flow. For convenience, we introduce some new notations 
G/yJRT 0 = (g x ,g y ), U/yJRT 0 = (u,v) and C/y/.RT 0 = c = (. c x ,c y ). This is just for the 
simplicity of formulations and should not be understood as a non-dimensional transformation 
in this work. Instead, we stick to the dimensional system presented in Eqs. (|T|) and Q. 
Moreover, for convenience, we use letter A to denotes the coefficient dt/{r + 0.5 dt) of 
the relaxation term and B for the coefficient rdt/{r + 0.5 dt) of the force term. Without 
influencing the discussion, we set g y and v to be zero. Therefore, for the D2Q4 lattice, the 
evolutionary rules for the jth bulk node are 


fj,i — (1 _ A)f j+ 14 + -Apuj + 1 + H- |A (13a) 

fj,2 = (1 - A)fj-1,2 + -^ApUj- 1 + H -jA (13b) 

Si,, = (1 - A)f j+ 1,3 - ±A n+1 + Y _ (13c) 

f'iA = (1 - A)li- 1,4 - \Apu^ + ^ - (13d) 

while Eq. ([ 6 ]) for the velocity becomes 


P u j ~ fj, 1 + fj ,2 “ fj ,3 - fj ,4 + 


PQxdt 


(14) 


To get the macroscopic governing equation, we follow the procedure of m For convenience, 


we write a few alternative variants of rules Eqs. (13)-(14), i.e., 


fj- 1,1 — (1 _ A)fjj + -Apuj + H-^A (15a) 

fj+ 1,2 = (1 ~ A)fj j2 + -Apuj + H-|A (15b) 

Si- 1,3 = (1 - + ^ - (15c) 

I j i M = (1 - 1 i/, :) - Fpip + 4^ - (15d) 
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and 


P u j+ 1 — fj+1,1 + fj+ 1,2 — /j+1,3 ~ /j+1,4 + 
PUj-i = fj- 1.1 + fj- 1,2 — fj—1,3 — fj- 1,4 + 


dtg x p 
2 ’ 


(16a) 

(16b) 


Applying the rules Eq. (13) into Eq. (14), we have 


prtj — (1 — A)(fj- 1,2 — fj- 1,4 + /j+i,i — fj+ 1,3) + + -Apuj + 1 + -BgxP + 2 ^' (-*- 1 ' 7 ) 

Hence, we need to work out (1 — A)(/j_i >2 — fj- 1,4 + /j+i,i — fj+ 1 , 3 ). To do so, the main 


idea is to use the rules Eqs.fll3h and (15) and Eqs. (14) and (16) alternatively. The aim is 


to relate the unknown distribution functions to the macroscopic quantities. For example, 


to obtain /j-1,2 — /j-1,4, we first use (16b) to represent the unknowns with Wj-i, /j- 1,1 and 


fj- i i3 , then use rules (15a) and (15c) to transform fj- 1,1 — fj- 1,3 into a formula of Uj , / J: 1 


and fj.'s- Finally, we can apply the rule Eq. (14) to convert all /)s into a formula of Uj. 
Following this idea and through a few iterations, we obtain 


Ag x p(Adt + 2 B) + (2 — A)pMj_! + 2(A — 2)p-Uj + (2 — A)puj + i = 0. 
Considering the meaning of A and B , the equation becomes 

pdt 2 g x + pr(uj_ 1 — 2?ij + ?i J+1 ) = 0. 

Further using dt = dx/^/RT 0 , p = pr, and po = P 0 RT 0 , the final form is 

p(Uj—\ — 2llj + Mj + i) 


dx 2 


+ P9x = 0, 


(18) 


(19) 


( 20 ) 


which is exactly the second central difference scheme of the NS equations for this simple 
force-driven flow. It has a simple solution 


(L-jdx)jdx . 

Uj = " - 2p -+ Ua ’ J = °> 2 ’ 3 ’ ’ ’ ’ N - 


( 21 ) 


The slip velocity is denoted by U s , which is produced by the boundary treatment, either 
physically or non-physically. So, its exact value will depend on the specific boundary condi¬ 
tion. 

To find U s , we also follow the procedure of [16]. First, we introduce a notation Uq = 
Yl a c a,xfo,a where the incoming distribution functions will be determined by the boundary 
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condition. In order to find the slip velocity, we actually look at the node j — 1. Following the 
manner of finding the bulk equation, we will be able to obtain the relation of the prescribed 
boundary speed u 0 , tq, u 2 and U 0 . The trick is that the prescribed boundary velocity (w 0 , 0) 
is used when applying Eq. (13) into Eq. 0- However, when using the rule ( 16b[ ) for j — 1, 
we need to consider the relation Uq = c a,xfo,a • Through simple calculations, the relation 
can be written as 

= up + u 2 + dt 2 g x _ (.dt - 2r) (Up - u 0 ) ^ 


2 ' 2 r 4 r 

For simplicity, we assume a non-slip boundary with zero speed at wall in the following, i.e., 


Uq is set to be zero. Then, substituting the solution (21) into Eq. (22), we can obtain the 
slip velocity 


U„ = 


2t — dt . 


2 T 


-Un. 


(23) 


For the so called modihed BB (MBB) scheme (means that collision and forcing still occur 
at boundary nodes [I2J), the rule at the boundary point is 


/o,4 — /o,i /o,2 — /o,3- 


(24) 


It can be seen that Uq must be zero. Hence, for the D2Q4 model, the slip velocity U s is zero. 
Although the SRT scheme is used, the MBB scheme leads to a correct non-slip boundary 
condition. Moreover, if rotating the wall direction from horizontal to vertical, it can be 
easily seen that slip velocity will also be zero. 

For the BB scheme without collision and forcing occurring at boundary nodes, we can 


not directly apply Eqs. (22) and (23). But it is straightforward to see that U s will be zero. 


As has been shown, the results of D2Q4 model are significantly different from that of 
the D2Q9 model and D2Q5 model presented in pT, 12]. This helps to clarify the relation 
between the SRT scheme and non-physical slip velocity. The SRT scheme plus the MBB or 
BB boundary scheme does not necessarily induces non-physical slip velocity. 


B. Mechanism of non-physical slip velocity with D2Q9 model. 

Now it is natural to ask why there is non-physical slip velocity in such as the D2Q9 
solution. For this purpose, we return to the D2Q9 model following |T2]- The rules for the 












D2Q9 model are 


f _ 4 P 
J j A g 


2 K 

9 


~ _ puj p Uj Bg x p p 
Ij > 2 9 3VS 3V3 A 9 

f j ,3 = ~YgAp u2 j- 1 + + (1 ~ A)fj- 1}3 

f = p A 

JjA g 


j- 

pUj _ Bg x p p 
3 a/3 3a/3A + 9 


fj ,5 - _ A).fj+\,5 


9 


A =h Apa J' + wr + *r + id + (1 - • 4)/j - 1 ’ 8 


f j J = l^ A P U U 


Apuj _ 

12\/3 


i-i Ap 
36 


-£>SLP 

12\/3 


; _ 1 „ 2 -4pwj+i , Ap 

h8 36 pUj+1 12 a/3 36 12 a/3 


+ (1 _ A )fj~ 1,7 

+ (1 _ A)fj+1,8 


f _ 1 A 2 , Apw j+ i , Ap , ^PxP , n ^ 

/ 3 ,„ - 3g^j + , + + 35 + ^ + (1 - W 


(25a) 

(25b) 

(25c) 

(25d) 

(25e) 

(25f) 

(25g) 

(25h) 

(25i) 


As discussed before, we will pay particular attention to discrete velocities parallel to the wall. 
With the horizontal wall, they are the 2nd and 4th velocity as shown in Fig. [l] Looking at 


Eqs. (25b) and (25d), we remind that they are mainly the consequence of periodic boundary 


conditions for the inlet and outlet, i.e., there is no gradient in the streamwise direction. In 
other words, they are not solely related to the SRT scheme and are NOT determined by the 
BB scheme at all. Similarly, the governing equation for bulk nodes is 


p(uj -1 — 2 Uj + u j+ 1) 
dx 2 


+ 3 dtlpg x - 0, 


(26) 


where dtg means the time step for the D2Q9 model. Assuming the space step is same for 
both two models, dtg = dt/\/ 3 where dt is time step for the D2Q4 model. Hence two models 
yield same governing equation for momentum. For brevity, we only discuss the MBB scheme. 
Therefore, U\ and U s can be written as 


u 0 + u 2 dt\g x 
U1 = —2~ + — 


3 (dtg - 2 t ) (Up - Up) 

At 


( 27 ) 


and 
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TT _ 3(2 r-dt 9 ) TT 

U S — -o- Uq. 

It 


( 28 ) 


It can be seen that the form of U s is consistent with Eq. (18) in |12| . To hnd Uq we need to 
calculate out 

/o ,2 — /o,4 + /o,6 — /o,7 + /o,9 ~ /o,8- (29) 

Applying the MBB rule 


/o,7 — /o,9 /o,6 — /o,8 /o,3 — /o,5, 


(30) 


we only need to consider / 0) 2 — fop- After simple calculations using Eqs. (25b) and (25d) 


(note that the MBB rule allows collisions at boundary, and again, these two equations 
are actually determined by the periodic boundary condition) with the prescribed boundary 
velocity (uq = 0, 0), we hnd it equals 


/o,2 — /o,4 — T^dxPT. 


Regarding that 


9x 


8pZl m 

Tv 


(31) 


(32) 


where u m denotes the centerline speed without slip velocity at boundaries, U s is written as 


U s = -g x (dt 9 - 2r) = 


8p(dtg — 2 r)u r 

Up 


(33) 


which is in the form of physical dimension. To transform to the commonly used lattice unit, 
we uses relations 


p = pRT 0 r L = Ndx = Ndtg \J 3 RT 0 


(34) 


and 


which yields 


r 1 
dtg + 2’ 


Us 


8(f - l)(2f - 1 )u m 
3N 2 


(35) 

(36) 


Here the units of U s and u m is not important since they cancel each other. We note the 


form Eq. (36) is slightly different from Eq. (22) in [12]. This is mainly because of difference 
of the factor B , i.e., the treatment of the body force term, cf. Eq. Q and Eq. (1) in [1(2 . 
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Clearly, the non-physical slip velocity obtained in |T 2 ] is due to the contribution of /o ,2 and 
/ 0 > 4 . However, as we have stressed, they are mainly the consequence of periodic boundaries 
at the streamwise direction. Therefore, the failure of the MBB scheme is due to lack of 
definition on the behaviour of discrete velocities parallel to wall. Then, they are actually 
controlled by other factors. In this case, it is the boundary scheme used in the inlet and 
outlet, which is not the bounce-back scheme. 

In this way, the slip velocity may be arbitrary in numerical practice, which may depend 
on the specific inlet and outlet boundary condition, geometry, other numerical operations 
at the boundary, and even the initial condition at wall. 


By identifying the mechanism of non-physical slip velocity, we may be able to devise 
remedy for the BB scheme. The key is to supplement the definition for the behaviour of the 
discrete velocities parallel to wall if there are any. Since other distribution function pairs 
can cancel each other when obtaining the velocity, they must also be able to cancel each 
other so that the velocity is zero. For instance, in this force-driven Poiseuille flow, although 
the bulk points must admit the consequence of the inlet and outlet boundary conditions, 
wall boundary points do not have to do so. In other words, we do not necessarily need to 


apply rules Eqs. (25b) and (25d) which has been done above and in m By contrary, We 
may initially set / 0) 2 and / 0i 4 to be a equilibrium distribution with zero velocity and they 
can remain their initial state all the time. This simple fix is able to correctly yield zero slip 
velocity. 


In practice, this can be incredibly easy for simple geometries. In the following section, 
we will show that actually no extra effort is necessary for a lid-driven cavity flow. However, 
the solution for complex geometries may need further investigation. 


On the other hand, it is worth noting that extra care may be necessary when using Eq. 


(36) to analyse the accuracy. At a first glance, U s seems to be a second order small quantity. 


However, in Eq.(33), it is actually more or less a constant error since we should set dtg to be 
smaller than r for stability while both g x and r are constant for a given incompressible and 
isothermal flow configuration. In our view, that is the confusion caused by using numerical 
time/spatial steps as reference quantities. 
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C. D2Q9 simulations for lid-driven cavity flow 


In this section, we will show how to utilise the above observation to devise remedy for the 
non-physical slip velocity. For this purpose, we will simulate the lid-driven cavity flow using 
the D2Q9 model. At the bottom, left and right wall, we will implement the bounce-back 
scheme. The non-equilibrium bounce-back scheme [18] is adopted for the top moving wall 
to bring in a wall velocity. 

For the cavity flow, we notice a fact that, for the discrete velocities parallel to wall, their 
distributions at boundary nodes are never affected by those of bulk. For the BB scheme, they 
will only be affected by their neighbours at wall. For the MBB scheme, the local collisions 
will also come into play. This fact can be utilised for eliminating the slip velocity. 

For the BB scheme, it can be easily seen that, for the discrete velocities parallel to wall 
(e.g., C 2 and C 4 for the bottom wall), the initial state is actually maintained in a way that 
the information is cycling among wall nodes. If the initial conditions at all wall nodes are 
set to be the uniform equilibrium distribution with zero velocity, the distribution of such as 
C -2 and C 4 will always be able to cancel each other when finding the velocity. Therefore, 
the non-slip velocity boundary can be achieved without any extra effort. The corner points 
at the top wall are a little more tricky as they are singular points. In practice, they may be 
treated as either a top wall point or a point of the left or right wall. Here, to maintain the 
benefit of “no extra” effort, we need to treat them as a left or right wall point. Otherwise, 
distributions at the left and right wall will be affected by those of the top wall which are 
changing with time, and the initial equilibrium state with zero velocity will break down. 
The other two corner points can treated in normal way although more discrete velocities 
need to be “bounced back”. 

For the MBB scheme, as collisions will occur at boundary points, the initial state cannot 
be maintained. However, using the fact that the information can not propagate into the bulk 
for the discrete velocities parallel to the wall, we are able to blend the relevant distribution 
to obtain the nonslip condition at wall. For instance, for the bottom wall, we can use the 
average of distribution of C 2 and C 4 as their new value after the streaming step. 

Numerical simulations are conducted for both two ways with four different Reynolds 
numbers while the top wall speed is fixed as 0.1^/ RTq. To examine the speed at the bottom, 
left, and right wall, we calculate the sum of \JU ■ U of all nodes at these three walls at every 
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Re = 10 

Re = 100 

Re = 500 

Re = 1000 

Average speed (no collision at wall, x 10 18 ) 

8.774 

8.740 

8.831 

8.844 

Average speed (with collision at wall, xlO -18 ) 

11.663 

11.594 

11.596 

11.610 


Table I. Average speed at the bottom, left, and right wall. The top wall speed is 0.1\/RTo. 


time step and then obtain the average speed per time step and per node. Since we are not 
examine the solution accuracy, no convergence test will be done for mesh size. By contrary, 
we will use as coarse mesh as possible to obtain results quickly. The maximum time step 
is set to be 10, 0000 iterations. While it is not of interest if the steady state is approached, 
the first order time derivative of the L 2 norm error of velocity is found to be smaller than 
3.5 x 10~ 4 except for cases with Re = 1000. The results are summarised in Table [1} As has 
been shown, the speed at walls are effectively zero within the machine resolution (double 
precision). It is worth noting again that actually no extra effort is necessary for the BB 
scheme. 

IV. CONCLUDING REMARKS 

To conclude, we have investigated the issue of the slip velocity at wall boundaries in 
lattice Boltzmann simulations with the BB scheme family. To identify the mechanism, we 
have analytically compared the solutions of a D2Q4 lattice and the commonly used D2Q9 
lattice for the force-driven Poiseuille flow. It is found that the BB family scheme does not 
define the behaviour of discrete velocities parallel to the wall. Mathematically the boundary 
condition is incompletely determined. This gives opportunities for other factors to affect 
them, such as the boundary conditions for inlet and outlet or even the initial condition at 
boundary. The non-physical slip velocity are induced exactly by these undesired effects. 
Therefore, the scheme for the bulk (e.g., the SRT scheme) is not the intrinsic reason for the 
slip velocity. 

To solve this issue, the key is to supplement the definition for the discrete velocities 
parallel to wall. By simulating the lid-driven cavity flows, We have shown that this can be 
incredibly easy for simple geometries. In fact, we may need no extra effort . The future 
study is to find if there is similar solution for complex geometries, which is already under 
progress. 
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